

In [2]:

    
from __future__ import division, print_function
import time
import os

# Third-party
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
from scipy.signal import argrelmin

# Custom
import gary.coordinates as gc
import gary.dynamics as gd
import gary.integrate as gi
import gary.io as io
import gary.potential as gp
from gary.units import galactic

from streammorphology.freqmap import estimate_dt_nsteps

Stream-fanning orbit



In [3]:

    
x0 = np.array([8.312877511, 0.242593717, 16.811943627])
v0 = ([-52.429087, -96.697363, -8.156130]*u.km/u.s).to(u.kpc/u.Myr).value
w0 = np.append(x0,v0)



In [5]:

    
t,w = potential.integrate_orbit(w0, dt=0.5, nsteps=12000, Integrator=gi.DOPRI853Integrator)



In [9]:

    
gd.peak_to_peak_period(t, w[:,0,2])









    Out[9]:





397.47954910714293



In [4]:

    
potential = gp.LM10Potential()



In [5]:

    
E = potential.total_energy(w0[:3], w0[3:])[0]

# find where it intersects X-Z plane
t,w = potential.integrate_orbit(w0, dt=0.5, nsteps=10000, Integrator=gi.DOPRI853Integrator)
ymin_ix = argrelmin(np.abs(w[:,0,1]))[0]
xz_ix = np.where((w[:,0,0] > 0) & (w[:,0,0] > 2))[0]

for i in ymin_ix:
    if i in xz_ix:
        print(i)
        break
        
print("Pal5 hits X-Z plane at: ({0:.2f},{1:.2f})".format(w[i,0,0], w[i,0,2]))
print("Energy: {0:.4f}".format(E))









    



5
Pal5 hits X-Z plane at: (8.18,16.79)
Energy: 0.0720



In [6]:

    
dt,nsteps = estimate_dt_nsteps(potential, w0, nperiods=5000, nsteps_per_period=512)
dt,nsteps









    Out[6]:





(0.4179181530510605, 5120000)



In [7]:

    
# path = "/Users/adrian/projects/morphology/output/pal5/"
# if not os.path.exists(os.path.join(path,"le.npy")):
#     a = time.time()
#     le,t,w = gd.fast_lyapunov_max(w0, potential, dt, nsteps)
#     print("Took {0:.2f} seconds".format(time.time() - a))

#     np.save(os.path.join(path,"le.npy"), le)
#     np.save(os.path.join(path,"t.npy"), t)
#     np.save(os.path.join(path,"w.npy"), w)

# le = np.load(os.path.join(path,"le.npy"))
# t = np.load(os.path.join(path,"t.npy"))
# w = np.load(os.path.join(path,"w.npy"))



In [8]:

    
a = time.time()
le,t,w = gd.fast_lyapunov_max(w0, potential, dt, nsteps, d0=1E-6)
print("Took {0:.2f} seconds".format(time.time() - a))









    



Took 3240.26 seconds



In [9]:

    
E = potential.total_energy(w[:,0,:3], w[:,0,3:])
plt.semilogy(np.abs((E[1:]-E[0])/E[0]), marker=None, drawstyle='steps')









    Out[9]:





[<matplotlib.lines.Line2D at 0x7f434f46d110>]



In [10]:

    
plt.figure(figsize=(10,8))
plt.loglog(t[1:-10:10], le, marker=None)









    Out[10]:





[<matplotlib.lines.Line2D at 0x7f4357058790>,
 <matplotlib.lines.Line2D at 0x7f4347d3ec10>]



In [22]:

    
# fig = gd.plot_orbits(w[::10,0], marker=',', alpha=0.5, linestyle='none')

Pal 5 Lyapunov time



In [25]:

    
1 / le[-1000:,0].mean()









    Out[25]:





4748.1778800226029

NAFF



In [22]:

    
dt,nsteps = estimate_dt_nsteps(potential, w0, nperiods=50, nsteps_per_period=512)
dt,nsteps









    Out[22]:





(0.4179181530510605, 50000)



In [23]:

    
t,w = potential.integrate_orbit(w0, dt=dt, nsteps=nsteps, Integrator=gi.DOPRI853Integrator)



In [60]:

    
r = np.sqrt(np.sum(w[:,0,:3]**2, axis=-1))
TT = gd.peak_to_peak_period(t, r)
12. / TT * 1000.









    Out[60]:





41.163572492775089



In [62]:

    
TT, T









    Out[62]:





(291.51988696088529, array([ 399.75006166,  427.85235185,  399.95648229]))



In [47]:

    
naff = gd.NAFF(t[:nsteps//2+1], p=2)



In [67]:

    
new_w = gc.cartesian_to_poincare_polar(w[:nsteps//2+1,0])
# new_w = w[:nsteps//2+1,0]
fs = [(new_w[:,i] + 1j*new_w[:,i+3]) for i in range(3)]
freqs1,d,nv = naff.find_fundamental_frequencies(fs)



In [68]:

    
new_w = gc.cartesian_to_poincare_polar(w[nsteps//2:,0])
# new_w = w[nsteps//2:,0]
fs = [(new_w[:,i] + 1j*new_w[:,i+3]) for i in range(3)]
freqs2,d,nv = naff.find_fundamental_frequencies(fs)



In [72]:

    
freqs1, freqs2









    Out[72]:





(array([-0.02157161,  0.01511077, -0.01570967]),
 array([-0.02155405, -0.01491793, -0.01569599]))



In [69]:

    
T = np.abs(2*np.pi / freqs1)
T









    Out[69]:





array([ 291.27109677,  415.80852345,  399.95648229])



In [70]:

    
nperiods = nsteps * dt / T
nperiods









    Out[70]:





array([ 71.74040914,  50.25367801,  52.24545314])



In [75]:

    
f = 0.005
R = np.abs((np.abs(freqs2)-np.abs(freqs1))/freqs1) / 25.
f / R









    Out[75]:





array([ 153.58279731,    9.79495191,  143.50893684])

KLD



In [25]:

    
ball_w0 = create_ball(w0, potential, N=1000, m_scale=2.5E4)
kld_dt,kld_nsteps = estimate_dt_nsteps(potential, w0, nperiods=50, nsteps_per_period=100)



In [26]:

    
kld_t,kld,mean_dens = do_the_kld(256, ball_w0, potential, kld_dt,kld_nsteps, 
                                 kde_bandwidth=10., density_thresholds=None)



In [29]:

    
# plt.loglog(kld_t / T.mean(), mean_dens)
plt.semilogx(kld_t, kld)
plt.axvline(6000)









    Out[29]:





<matplotlib.lines.Line2D at 0x1105d5890>



In [30]:

    
plt.loglog(kld_t, mean_dens)
plt.axvline(6000)









    Out[30]:





<matplotlib.lines.Line2D at 0x110a2de90>



In [ ]: